#include "channels.hpp"

channels::~channels(){}

channels::channels(){
  dim_of_H=0;
  number_of_channels=0;
}

channels::channels(const channels& chan){
  if(this != &chan){
    dim_of_H=chan.dim_of_H;
    number_of_channels=chan.number_of_channels;
    ops=chan.ops;
    D_rates=chan.D_rates;
    L_rates=chan.L_rates;
  }
}

channels& channels::operator=(const channels& chan){
  if(this != &chan){
    dim_of_H=chan.dim_of_H;
    number_of_channels=chan.number_of_channels;
    ops=chan.ops;
    D_rates=chan.D_rates;
    L_rates=chan.L_rates;
  }
  return *this;
}



void channels::set_op(matrix<double> o, int ii){
  if(dim_of_H == 0){
    dim_of_H=o.size1();
  }
  //If channel with index ii exists
  if( ops.find(ii) != ops.end()){
    ops[ii]=o;
  }
  else{
    ops[ii]=o;
  }
}

void channels::set_op(matrix<double> o){
  if(dim_of_H == 0){
    dim_of_H=o.size1();
  }
  
  //Key of the one past last value
  int last=ops.end()->first;
  ops[last]=o;
}

void channels::set_decay_rate(rate r, int ii){
  //If channel with index ii exists
  if( D_rates.find(ii) != D_rates.end()){
    D_rates.erase(ii);
    D_rates.insert(std::make_pair(ii,r));
  }
  else{
    D_rates.insert(std::make_pair(ii,r));
  }
}

void channels::set_decay_rate(rate r){
  int last=D_rates.end()->first;
  D_rates.insert(std::make_pair(last,r));
}

void channels::set_lamb_rate(rate r, int ii){
  //If channel with index ii exists
  if( L_rates.find(ii) != L_rates.end()){
    L_rates.erase(ii);
    L_rates.insert(std::make_pair(ii,r));
  }
  else{
    L_rates.insert(std::make_pair(ii,r));
  }
}

void channels::set_lamb_rate(rate r){
  int last=L_rates.end()->first;
  L_rates.insert(std::make_pair(last,r));
}

int channels::get_number_of_channels(){
  int n_ops=ops.size();
  int n_rates=D_rates.size();
  if(n_ops>n_rates){
    std::cout<<"Channels not correctly specified"<<std::endl;
    std::cout<<"Number of operators "<< n_ops
	     <<"is larger than number of rates "<<n_rates<<"."<<std::endl;
    exit(1);
  }
  if(n_rates>n_ops){
    std::cout<<"Channels not correctly specified"<<std::endl;
    std::cout<<"Number of operators "<< n_ops
	     <<"is smaller than number of rates "<<n_rates<<"."<<std::endl;
    exit(1);
  }
  return n_ops;
}
 
rate channels::get_decay_rate(int ii){
  return D_rates[ii];
}

rate channels::get_lamb_rate(int ii){
  return L_rates[ii];
}
 
matrix<double> channels::get_op(int ii){
  return ops[ii];
}
matrix<double> channels::operator()(int ii,double t){
  return D_rates[ii](t)*ops[ii];
}
matrix<double> channels::operator()(double t){
  matrix<double> mat(dim_of_H,dim_of_H,0);
  if( !L_rates.empty() ){
    if(ops.size()==L_rates.size()){
      for(std::pair< std::map<int,matrix<double> >::iterator,std::map<int,rate>::iterator >
	    ii(ops.begin(),L_rates.begin());
	  ii.first != ops.end() && ii.second != L_rates.end();
	  ++ii.first,++ii.second){
	double r=(*ii.second).second(t);
	matrix<double> m=(*ii.first).second;
	mat += r*prod(herm(m),m);
	
      }
    }
    else{
      std::cout<<"Number of Lamb rates is not zero and \n"
	       <<"it does not match with the number of jump operators."
	       <<std::endl;
      exit(1);
    }
  }
  return mat;
}


matrix<double> channels::get_lamb_shift(double t){
  return this->operator()(t);
}
